* Creates Figure 2

clear all

use "$savedata/masterdata.dta", replace

keep if sample25==1
gen vol = vol25

preserve

foreach x in 30{
* Residuals with no controls
reghdfe survive`x', absorb(i.hyid i.dow##i.admidate_mont##i.finyear, savefe) keepsingleton
predict residual, residuals

xtreg residual, fe i(doctor_id)

gen sigma_u = e(sigma_u)
gen sigma_e = e(sigma_e)
predict docfe`x', u

gen signal_2step = ((sigma_u^2) / ((sigma_u^2) + ((sigma_e^2)/vol)))
gen adj_docfe`x' = docfe`x'*signal_2step

* Residuals with all controls
drop sigma* residual* signal*
reghdfe survive`x' c.prevyear_cost c.std_ami3 c.std_nonami3, absorb(i.sex##i.derv_age i.black i.mixed i.chinese i.asian i.race_miss i.ynch* i.di1 i.di2 i.di3 i.di4 i.di5 i.shock i.arythmia i.arthero i.arrest i.hyid i.dow##i.admidate_mont##i.finyear, savefe) keepsingleton
predict residual, residuals

xtreg residual, fe i(doctor_id)
predict docfe`x'_control, u
gen sigma_u = e(sigma_u)
gen sigma_e = e(sigma_e)

gen signal_2step = ((sigma_u^2) / ((sigma_u^2) + ((sigma_e^2)/vol)))
gen adj_docfe`x'_control = docfe`x'_control*signal_2step


* Predicted residuals - first need to predict mortality, and then use as outcome

reghdfe survive`x' c.prevyear_cost c.lsoa_imdscore10 c.lsoa_wa_benefits c.sales c.prevyear_admit_em c.prevyear_admit_el, absorb(i.derv_age##i.sex i.black i.mixed i.chinese i.asian i.race_miss i.ynch* i.prevyear_stroke i.di1 i.di2 i.di3 i.di4 i.di5 i.shock i.arrest i.arythmia i.arthero, savefe) keepsingleton 
predict pred_survive`x', xbd
egen std_pred_survive`x' = std(pred_survive`x')

** Now use predicted value to get residuals

drop sigma* residual* signal*
* Residuals with all controls
reghdfe pred_survive`x', absorb(i.hyid i.dow##i.admidate_mont##i.finyear, savefe) keepsingleton
predict residual, residuals

xtreg residual, fe i(doctor_id)

gen sigma_u = e(sigma_u)
gen sigma_e = e(sigma_e)
predict pred_docfe`x', u

gen signal_2step = ((sigma_u^2) / ((sigma_u^2) + ((sigma_e^2)/vol)))
gen adj_pred_docfe`x' = pred_docfe`x'*signal_2step

egen std_adj_pred_docfe`x' = std(adj_pred_docfe`x')

drop sigma* residual* signal*
}



collapse (mean) adj* (sum) one, by(pconsult)

foreach x in 30{
xtile n20_`x' = adj_docfe`x', n(20)
xtile n20_p_`x' = adj_pred_docfe`x', n(20)
xtile n20_c_`x' = adj_docfe`x'_c, n(20)

}

collapse (mean) adj* (mean) one, by(n20_c_30)

forval z=1(1)20{
	
	su adj_pred_docfe30 if n20_c_30==`z'
	matrix observe_pred`z' = r(mean)
		
}

forval z=1(1)20{
	
	su adj_docfe30 if n20_c_30==`z'
	matrix observe_noc`z' = r(mean)
		
}

restore

capture program drop myboot2

program define myboot2, rclass

* draw my bootstrap sample
preserve
bsample, strata(doctor_id)

foreach x in 30{
* Residuals with no controls
reghdfe survive`x', absorb(i.hyid i.dow##i.admidate_mont##i.finyear, savefe) keepsingleton
predict residual, residuals

xtreg residual, fe i(doctor_id)

gen sigma_u = e(sigma_u)
gen sigma_e = e(sigma_e)
predict docfe`x', u

gen signal_2step = ((sigma_u^2) / ((sigma_u^2) + ((sigma_e^2)/vol)))
gen adj_docfe`x' = docfe`x'*signal_2step

* Residuals with all controls
drop sigma* residual* signal*
reghdfe survive`x' c.prevyear_cost c.std_ami3 c.std_nonami3, absorb(i.sex##i.derv_age i.black i.mixed i.chinese i.asian i.race_miss i.ynch* i.di1 i.di2 i.di3 i.di4 i.di5 i.shock i.arythmia i.arthero i.arrest i.hyid i.dow##i.admidate_mont##i.finyear, savefe) keepsingleton
predict residual, residuals

xtreg residual, fe i(doctor_id)
predict docfe`x'_control, u
gen sigma_u = e(sigma_u)
gen sigma_e = e(sigma_e)

gen signal_2step = ((sigma_u^2) / ((sigma_u^2) + ((sigma_e^2)/vol)))
gen adj_docfe`x'_control = docfe`x'_control*signal_2step

* Predicted residuals - first need to predict mortality, and then use as outcome

reghdfe survive`x' c.prevyear_cost c.lsoa_imdscore10 c.lsoa_wa_benefits c.sales c.prevyear_admit_em c.prevyear_admit_el, absorb(i.derv_age##i.sex i.black i.mixed i.chinese i.asian i.race_miss i.ynch* i.prevyear_stroke i.di1 i.di2 i.di3 i.di4 i.di5 i.shock i.arrest i.arythmia i.arthero, savefe) keepsingleton 
predict pred_survive`x', xbd
egen std_pred_survive`x' = std(pred_survive`x')

** Now use predicted value to get residuals

drop sigma* residual* signal*
* Residuals with all controls
reghdfe pred_survive`x', absorb(i.hyid i.dow##i.admidate_mont##i.finyear, savefe) keepsingleton
predict residual, residuals

xtreg residual, fe i(doctor_id)

gen sigma_u = e(sigma_u)
gen sigma_e = e(sigma_e)
predict pred_docfe`x', u

gen signal_2step = ((sigma_u^2) / ((sigma_u^2) + ((sigma_e^2)/vol)))
gen adj_pred_docfe`x' = pred_docfe`x'*signal_2step

egen std_adj_pred_docfe`x' = std(adj_pred_docfe`x')

drop sigma* residual* signal*
}

collapse (mean) adj* (sum) one, by(pconsult)

foreach x in 30{
xtile n20_`x' = adj_docfe`x', n(20)
xtile n20_p_`x' = adj_pred_docfe`x', n(20)
xtile n20_c_`x' = adj_docfe`x'_c, n(20)

}

collapse (mean) adj* (mean) one, by(n20_c_30)

forval z=1(1)20{
	
	su adj_pred_docfe30 if n20_c_30==`z'
	return scalar d_pred`z' = r(mean)
		
}

forval z=1(1)20{
	
	su adj_docfe30 if n20_c_30==`z'
	return scalar d_noc`z' = r(mean)
		
}


restore

end


** Then run bootstrap here

#delimit ;
* Simulate 199 times to test;
simulate diff_pred1 = r(d_pred1) diff_pred2 = r(d_pred2) diff_pred3 = r(d_pred3) diff_pred4 = r(d_pred4) diff_pred5 = r(d_pred5)
diff_pred6 = r(d_pred6) diff_pred7 = r(d_pred7) diff_pred8 = r(d_pred8) diff_pred9 = r(d_pred9) diff_pred10 = r(d_pred10)
diff_pred11 = r(d_pred11) diff_pred12 = r(d_pred12) diff_pred13 = r(d_pred13) diff_pred14 = r(d_pred14) diff_pred15 = r(d_pred15)
diff_pred16 = r(d_pred16) diff_pred17 = r(d_pred17) diff_pred18 = r(d_pred18) diff_pred19 = r(d_pred19) diff_pred20 = r(d_pred20)

diff_noc1 = r(d_noc1) diff_noc2 = r(d_noc2) diff_noc3 = r(d_noc3) diff_noc4 = r(d_noc4) diff_noc5 = r(d_noc5)
diff_noc6 = r(d_noc6) diff_noc7 = r(d_noc7) diff_noc8 = r(d_noc8) diff_noc9 = r(d_noc9) diff_noc10 = r(d_noc10)
diff_noc11 = r(d_noc11) diff_noc12 = r(d_noc12) diff_noc13 = r(d_noc13) diff_noc14 = r(d_noc14) diff_noc15 = r(d_noc15)
diff_noc16 = r(d_noc16) diff_noc17 = r(d_noc17) diff_noc18 = r(d_noc18) diff_noc19 = r(d_noc19) diff_noc20 = r(d_noc20)

, reps(199) seed(32786105): myboot2;


#delimit ;
* Boostrap;
bstat, stat(observe_pred1, observe_pred2, observe_pred3, observe_pred4, observe_pred5, observe_pred6, observe_pred7, observe_pred8, observe_pred9, observe_pred10, 
observe_pred11, observe_pred12, observe_pred13, observe_pred14, observe_pred15, observe_pred16, observe_pred17, observe_pred18, observe_pred19, observe_pred20, 

observe_noc1, observe_noc2, observe_noc3, observe_noc4, observe_noc5, observe_noc6, observe_noc7, observe_noc8, observe_noc9, observe_noc10, 
observe_noc11, observe_noc12, observe_noc13, observe_noc14, observe_noc15, observe_noc16, observe_noc17, observe_noc18, observe_noc19, observe_noc20);

#delimit cr
* Output results
putexcel set "$results/vintiles_ses.xlsx", replace
putexcel A1 = "Variable"
putexcel B1 = "Coefficient"
putexcel C1 = "Std. error"

matrix b = e(b)'
putexcel A2 = matrix(b), rownames

matrix se = e(se)'
putexcel C2 = matrix(se)


* make into a figure

svmat b, name(coef)
svmat se, name(se)

keep coef se

gen n = _n if _n<21

gen coef2 = .
gen se2 = .
forval z=1(1)20{
	
	replace coef2 = coef1[_n+20] if _n==`z'	
	replace se2 = se1[_n+20] if _n==`z'	
}

drop if n>20

foreach x in 1 2{
gen lb`x' = coef`x' - 1.96*se`x'
gen ub`x' = coef`x' + 1.96*se`x'
}

foreach var in coef1 coef2 lb1 lb2 ub1 ub2{
	replace `var'=`var'*100
}

save "$savedata/bootstrap_s30.dta", replace

twoway scatter coef1 coef2 n, msymbol(dot square) mcolor(black red) msize(small small) || rcap lb1 ub1 n, lcolor(black) || lpoly coef1 n, lcolor(black) || rcap lb2 ub2 n, lcolor(red) || lpoly coef2 n, lcolor(red) ytitle("Percentage points") xtitle("Ventiles of doctor quality")  legend(order(1 "Predicted Survival" 2 "Observed Survival (no controls)") ) graphregion(color(white)) yscale(range(-4(1)4)) ylabel(-4(1)4)
graph export "$results/fig2.pdf", as(pdf) replace
